This code chunk takes a long time to compute but needs to be run once to generate the .RData file with individual data fits. The file was larger than 100mb so we couldn’t upload it to github.
# participant 45 is removed from this analysis because they gave a 0 surprise rating
# across all trials
df.fit = df.exp4.regression %>%
filter(participant != 45)
# restrict priors to be positive
priors = prior(normal(0, 10), class = b, lb = 0)
# initial model fits (used for compilation)
fit.brm_exp4_responsibility_baseline_individual =
brm(formula = responsibility ~ 1,
data = df.fit %>%
filter(participant == 1),
cores = 4,
chains = 4,
iter = 4000,
seed = 1,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = T),
file = str_c("cache/fit_brm_exp4_responsibility_baseline_individual"))
fit.brm_exp4_responsibility_importance_individual =
brm(formula = responsibility ~ 1 + importance,
data = df.fit %>%
filter(participant == 1),
prior = priors,
cores = 4,
chains = 4,
iter = 4000,
seed = 1,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = T),
file = str_c("cache/fit_brm_exp4_responsibility_importance_individual"))
fit.brm_exp4_responsibility_surprise_individual =
brm(formula = responsibility ~ 1 + surprise,
data = df.fit %>%
filter(participant == 1),
prior = priors,
cores = 4,
chains = 4,
iter = 4000,
seed = 1,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = T),
file = str_c("cache/fit_brm_exp4_responsibility_surprise_individual"))
fit.brm_exp4_responsibility_importance_surprise_individual =
brm(formula = responsibility ~ 1 + importance + surprise,
data = df.fit %>%
filter(participant == 1),
prior = priors,
cores = 4,
chains = 4,
iter = 4000,
seed = 1,
control = list(adapt_delta = 0.95),
save_pars = save_pars(all = T),
file = str_c("cache/fit_brm_exp4_responsibility_importance_surprise_individual"))
# update model fits for all participants
df.exp4.responsibility_individual_fit = df.fit %>%
group_by(participant) %>%
nest() %>%
ungroup() %>%
# fit model for each participant
mutate(fit_baseline = map(
.x = data,
.f = ~ update(fit.brm_exp4_responsibility_baseline_individual,
newdata = .x,
seed = 1)),
fit_importance = map(
.x = data,
.f = ~ update(fit.brm_exp4_responsibility_importance_individual,
newdata = .x,
seed = 1)),
fit_surprise = map(
.x = data,
.f = ~ update(fit.brm_exp4_responsibility_surprise_individual,
newdata = .x,
seed = 1)),
fit_importance_surprise = map(
.x = data,
.f = ~ update(fit.brm_exp4_responsibility_importance_surprise_individual,
newdata = .x,
seed = 1))) %>%
mutate(fit_baseline = map(fit_baseline, ~ add_criterion(.,
criterion = "loo",
moment_match = T)),
fit_importance = map(fit_importance, ~ add_criterion(.,
criterion = "loo",
moment_match = T)),
fit_surprise = map(fit_surprise, ~ add_criterion(.,
criterion = "loo",
moment_match = T)),
fit_importance_surprise = map(fit_importance_surprise,
~ add_criterion(.,
criterion = "loo",
moment_match = T)),
r_importance = map2_dbl(.x = data,
.y = fit_importance,
.f = ~ cor(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
r_surprise = map2_dbl(.x = data,
.y = fit_surprise,
.f = ~ cor(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
r_importance_surprise = map2_dbl(.x = data,
.y = fit_importance_surprise,
.f = ~ cor(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
rmse_baseline = map2_dbl(.x = data,
.y = fit_baseline,
.f = ~ rmse(.x$wrongfulness, .y %>%
fitted() %>%
.[, 1])),
rmse_importance = map2_dbl(.x = data,
.y = fit_importance,
.f = ~ rmse(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
rmse_surprise = map2_dbl(.x = data,
.y = fit_surprise,
.f = ~ rmse(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
rmse_importance_surprise = map2_dbl(.x = data,
.y = fit_importance_surprise,
.f = ~ rmse(.x$responsibility, .y %>%
fitted() %>%
.[, 1])),
model_comparison = pmap(.l = list(baseline = fit_baseline,
importance = fit_importance,
surprise = fit_surprise,
importance_surprise = fit_importance_surprise),
.f = ~ loo_compare(..1, ..2, ..3, ..4)),
best_model = map_chr(.x = model_comparison,
.f = ~ rownames(.) %>% .[1]),
best_model = factor(best_model,
levels = c("..1", "..2", "..3", "..4"),
labels = c("baseline",
"importance",
"surprise",
"importance_surprise")))
save(list = c("df.exp4.responsibility_individual_fit"),
file = "data/exp4_responsibility_individual_fit.RData")